Kulminacją dyskusji na temat tempa reakcji jądrowych jest wstawienie wszystkich wcześniejszych wyników: współczynnika tunelowania oraz S-factora do całki uśredniającej względem rozkładu Maxwella w temperaturze T.  Referencyjny wynik w postaci funkcji interpolacyjnej (fitu) to:

pp=4.01*10^-15 T9^(-2/3) Exp[-3.380/T9^(1/3)]*(1+0.123 T9^(1/3)+1.09*T9^(2/3)+0.938 T9)

S[En_]:=S0+Sprim*En

$Assumptions=kT>0&&Egamov>0

Całka jest możliwa do wyznaczenia analitycznie przy użyciu funkcji MeijerG (UWAGA: Mathematica czasem liczy tą funkcję, jej wartości numeryczne,  ŹLE ! ).

σv=Sqrt[8/(Pi (mp/2) (kT)^3)]Integrate[S[En]Exp[-Sqrt[Egamov/En]] Exp[-En/kT],{En,0,Infinity},Assumptions->kT>0]

Natomiast uświęcony przez czas sposób jej liczenia polega na zastosowaniu metody Laplace’a, czyli zamianie funkcji w kształcie piku na f. Gaussa. Szukamy maksimum:

E0=En/.Solve[D[Exp[-Sqrt[Egamov/En]] Exp[-En/kT],En]==0,En][[1]]

Przybliżenie kwadratowe na szczycie piku:

parabola=Series[-Sqrt[Egamov/En]-En/kT,{En,E0,2}]//Normal//PowerExpand

max=parabola/.En->E0//Exp

Jakość przybliżenia nie jest imponująca, ale ma ono działać dla małych temperatur. Widoczny na wykresie pik, nazywamy pikiem Gamowa.

Block[{Egamov=1,kT=Egamov/30},Plot[{Exp[-Sqrt[Egamov/En]] Exp[-En/kT],Exp[parabola]},{En,0, 1/2Egamov},PlotRange->{0,max}]]

Optymalna energia cząstki (protonu) przy której zachodzi najwięcej reakcji jest wynikiem 2 konkurujących czynników: 
1) rosnącego wraz z energią współczynnika tunelowania,
2) malejąca zgodnie z rozkładem Boltzmanna  liczba cząstek o danej energii.

Block[{Egamov=1,kT=Egamov/30},Plot[{Exp[-Sqrt[Egamov/En]] ,Exp[-En/kT]},{En,0, 1/2Egamov}]]

Rozszerzamy przedział całkowania na całą oś liczbową i voila!

σvSaddle=Sqrt[8/(Pi (mp/2) (kT)^3)]Integrate[S[En]Exp[parabola],{En,-Infinity,Infinity}]

σvSaddle=Sqrt[8/(Pi (mp/2) (kT)^3)]Exp[-((Egamov^(1/3))/((2^(2/3)) (kT^(1/3))))-((2^(1/3)) (Egamov^(1/3)))/(kT^(1/3))]Integrate[S[En+E0]Exp[-((3 (((En))^2))/(2 (2^(1/3)) (Egamov^(1/3)) (kT^(5/3))))],{En,-Infinity,Infinity}]

%/%%//Simplify

ε0=8.8541878128*10^-12;

barn=10^-28; (* m^2 *)

e=(801088317/5000000000000000000000000000);

keV=e*1000;

S0=4.07*10^-22*keV*barn

Sprim=4.52*10^-24*barn

kB=(1380649/100000000000000000000000000000);

h=(132521403/200000000000000000000000000000000000000000);

mp=1.6726219114600650621071552390089`10.*^-27;

NA=602214076000000000000000;

c=299792458;

Egamov=((4 (π^2) (e^4) mp)/(16 ((ε0)^2) (h^2)))

Egamov2=mp c^2*(Pi/137.036)^2 (* https://en.wikipedia.org/wiki/Gamow_factor *)

{Egamov2,Egamov}/keV

λpp=Simplify[σv/.{kT->kB*T9*10^9},T9>0]

λppSaddle=Simplify[σvSaddle/.{kT->kB*T9*10^9},T9>0]

pp//Simplify

{λpp,λppSaddle}/.T9->0.015

LogLogPlot[{λpp,λppSaddle,pp/NA*10^-6},{T9,0.01,10},PlotStyle->{Red,Blue,Black}]

{λpp,pp/NA}/.T9->0.015

Block[{T9=0.015,kT=kB*T9*10^9},Sqrt[8/(Pi (mp/2) (kT)^3)]NIntegrate[S[En]Exp[-Sqrt[Egamov/En]] Exp[-En/kT],{En,0,Egamov,Infinity}]]

Block[{ζ=1,ρ=1},Table[Exp[0.188 ζ Sqrt[ρ] (T9/1000)^(-3/2)],{T9,0.01,0.1}]]